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Abstract 

We have implemented the linear response approximation of a method proposed to compute the 
electron transport through correlated molecules based on the time-independent Wigner function 
[P. Delaney and J. C. Greer, Phys. Rev. Lett. 93, 036805 (2004)]. The results thus obtained for the 
zero-bias conductance through a quantum dot both without and with correlations demonstrate that 
this method is neither quantitatively nor qualitatively able to provide a correct physical description 
of the electric transport through nanosystems. We present an analysis indicating that the failure 
is due to the manner of imposing the boundary conditions, and that it cannot be simply remedied. 



PACS numbers: 73.63.-b, 73.63.Kv, 73.23.-b, 73.23.Hk, 71.10.Fd 



I. INTRODUCTION 



Electronic transport in artificial nanosystems and single molecules represents a topic of 
continuing current interest and remains a challenge both for fundamental science and tech- 
nological applications. In spite of numerous advances, there still exist important issues 
which are far from being resolved. Notoriously, the comparison between experimental and 
theoretical values for the current flowing through molecules usually show large discrepancies, 
typically several orders of magnitudes.-^ No consensus has been reached so far whether 
the discrepancies are to be attributed to uncontrollable experimental factors (compare mea- 
surements on same systems in Refs. Q^lf]!?]) or inadequate theoretical frameworks. 

Early theoretical calculations were carried out within the Landauer formalism,-^ which is 
based on a single-particle description. Because systematically the values of the currents thus 
obtained are at odds with those experimentally measured, a series of theoretical methods 
was developed to treat electron correlations, which are excluded by the Landauer approach. 

The description of electron correlation represents an important challenge from the theoret- 
ical side. The most popular theoretical approaches for nonequilibrium transport in correlated 
nanoscopic/molecular systems are based on nonequilibrium Green functions (NEGF),— ^ 
time-dependent DMRG (density matrix renormalization group),— and numerical renor- 
malization group (NRG).— In ab initio modeling of molecular electronics, the NEGF tech- 
nique is by far the most utilized one. It is usually combined with electron structure calcula- 
tions based on density functional theories (DFT).— The NEGF + DFT approach is con- 
ceptually simple and computationally less demanding than many-body methods.— i22i2L22i22 
However, because of the incompletely elucidated current dependences^ and self-interaction 
corrections^^ of the local exchange and correlation functionals, the DFT currents can sen- 
sitively vary by choosing different approximate potentials. 
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seems to be one of the 



Out of the methods proposed so far, the method proposed in Ref. 
most appealing, because of its claim to correctly reproduce the steady-state current experi- 
mentally measured through correlated molecular systems. It relies upon genuine many-body 
calculations and is not limited to a particular configuration interaction (CI), e. g., restricted 
to singly, doubly, or triply-excited configurations. The originality of this approach consists 
in the manner of imposing the boundary conditions. Rather than the Fermi distribution 
function, the Wigner distribution function is employed instead for the constrains at the 
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device-electrodes (reservoirs) interfaces. Unlike the former, which is meaningful only within 
a single-particle picture, the latter can be directly used in a many-body approach. Because 
the key quantity of this method is the stationary Wigner function (SWF), hereafter it will 
be referred to as the SWF-method. 

One should emphasize from the very beginning that the support of the SWF-method, and 
especially its manner to impose boundary conditions, entirely relies upon its ability to pro- 
vide currents comparable with those measured experimentally for two molecular systems^^ 
and plausible results for another class of molecules.— Because this method is conceptually so 
different from the other widely utilized approaches, it would be highly desirable to inquire 
the validity of this method within the realm of theory. If and to the extent to which the 
SFW-method is able to reproduce well established results, one can consider its predictions 
reliable. It is the purpose of the present work to investigate whether the SWF-method is 
able to correctly describe nanosystems whose behavior is well understood. 

The remaining part of the paper is organized in the following manner. The theoretical 
framework, the linear response approximation of the SWF-method, will be developed in Sec. 
HTl The model for the "device" we shall consider, consisting of a single quantum dot (QD), 
will be exposed in Sec. IIHI Next, in Sec. IIVt we shall compare the zero-bias conductance 
calculated without correlations by means of the SWF-method with the exact one. Sec. |V] 
will be devoted to the linear transport through the QD in the presence of correlations. The 
implications of our results will be discussed in Sec. ED Conclusions will be presented in Sec. 



II. METHOD 



Theoretical studies of the transport in nanosystems are inherently faced with the problem 
of separating the total system (S) into left and right reservoirs (electrodes), and device. The 
contacts (interfaces) between device and reservoirs (located at qL,R) are necessarily subject 
to arbitrariness. Depending on how demanding the numerical calculations are, parts of the 
electrodes (as large as possible) are included in the central part and treated as accurate as 
possible. In the SFW-method, the challenge is to determine the wave function \l/ for the 
total system S. 



Following Ref. |20|, we shall consider the many-electron wave function \l/ for the total 
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system S and its associate exact energy 8 in the presence of an applied bias 



£ = {V\H T \V) = (V\H\V) + W\V), (1) 

where H is the Hamiltonian of the total system S, and W is the perturbation due to applied 
bias. 

Within the single-particle picture, in a transport problem the (semi-infinite) reservoirs 
fix the electron chemical potentials \il,r at the left and right contacts, and the current flow 
is driven by the imbalance fiL — = created by an applied voltage. At either contact, 
the electron distribution is dictated by the reservoir Fermi function with the correponding 
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is that 



chemical potential. The attractive feature of the approach proposed in Ref. 
it is based on the Wigner function, which is directly applicable to a many-body system, 
without the need to resort to single-particle Fermi distribution functions. For problems 
describable within a single-particle picture, the Wigner function method to account for 
open boundaries specific for transport problems has certain advantages over the conventional 
scattering approach.— However, these advantages are merely technical there. 

For a system characterized by a many-electron wave function the one-particle Wigner 
distribution function f(q,p) is defined by 

f(q,p) = £ <tf| al^ a a q+r , a |tf > e'^ n . (2) 

r,<r 

Unlike the methods based on the NEGF, where semi-infinite leads are accounted for, the 
SWF-method was devised to be suited for ab initio quantum chemical calculations, where 
rather than semi-infinite electrodes, only (usually very) small parts thereof can be dealt with. 
Instead of having an infinite extension, the wave function \l/ and the related summation over 
r of Eq. (T3J) are inherently limited in space. In this paper, we shall adopt a one-dimensional 
discrete representation, wherein in a lattice of size N = 2M + 1 the site index ranges from 
— M to +M. In Eq. (T5]), a^ a (a' lcr ) denote the creation (annihilation) operators for electrons 
of spin a at site I. 

The open boundary conditions to be imposed should distinguish between electrons flowing 
into and those flowing out of the device.— The distribution of the former should be deter- 
mined by the reservoirs. In terms of the Wigner function, this means that, at the left (qi,) 
and right (q R ) boundaries between device and reservoirs, / (qL,PL > 0) and / (qn, Pr < 0) 
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should be fixed at values dictated by the reservoirs 3 



/ (Ql,Pl) = Sl (Pl, Hl) , for p L > 0, (3) 
/ (Qr, Pr) = fit (Pr; Vr) , for PR < 0. (4) 

These constraints are physically plausible and have been successfully applied previously to 
problems treated within single particle approximations.— i 3 - 1 -* 3 ^ The r.h.s. of Eqs. (J3j) and (J4]) 
represent the Fermi functions of the reservoirs with the chemical potentials shifted by the 
bias voltage, Hl,r — n±eV/2. Because this procedure cannot be applied for the many-body 
case, to determine the above Jl,r, which are basically properties of the reservoirs, it was 
claimed^ that they can be extracted from the wave function of the total system S at zero 
bias. That is, at zero temperature, one has to determine the ground state \&o of H, then 
calculate f by using \l/ instead of \I/ in Eq. (j2J), and impose 

/ (Ql,Pl) = fo (Ql,Pl) , for Pl > 0, (5) 
/ (Qr,Pr) = fo (Qr-i Pr) , for Pr < 0. (6) 

To describe a steady state, one has to impose in addition the condition that the average 
J of the electric current operator (e is the elementary charge) 33 

* = it ^H { a l+i,v a i,° ~ a {* a i+i,o) ( 7 ) 

be site (I) independent 

J, = (*|j,|*) = J. (8) 
Above, ti denotes the hopping integral (— M < I < M — 1). 



According to Ref. |2Cj, the solution of the transport problem is obtained by minimizing £ 



with the supplementary constraints (JSJ), (J6]), and © along with the normalization condition 

(* |^) = 1. (9) 

In this way, to reach its steady state, the device is allowed to optimize the Wigner distribution 
function of the electrons inside the device and that of the electrons flowing from it into 
reservoirs. (The distribution of outgoing electrons should depend only of the state of the 
device.) 

In the present Section, we shall work out the linear response approximation of the SFW- 
method, which enables us later to compute the zero-bias conductance. That is, we shall only 



consider changes to the relevant quantities of the order 0(V), caused by a small applied volt- 
age V. The system (left reservoir, device, and right reservoir), which consists of a collection 
of discrete sites ranging from — M to +M, is perturbed by a small electric perturbation 

M 

W = -e Yl nm, (10) 

l=-M 

where rii = a \a a i,° ls the electron number operator and <pi is the electric potential at site I. 
{f-M — <f+M — ~V/2). Starting with a system in the ground state \l/o; we shall look 

for a solution \l/ describing a steady state that slightly differs from \l/o- The wave function 
\l/ will be expanded in terms of the complete set of eigenstates of H (H\^f n ) = E n \^ n )) 

|*)=4>|*0> + £4»|*n>, (11) 

where A n = 0(V) for n ^ 0. One should note at this point that, although the eigenstates 
ty n can and will be chosen real, in order to satisfy the boundary conditions ([5} and ([6]) the 
expansion coefficients A n must be complex. While in general the minimization of 8, Eq. (JT]), 
with the constraints (jSJ), (EI), (jHJ), and (EJ) represents a nonlinear problem, it considerably 
simplifies in the linear response limit. By introducing the Lagrange multipliers Xi,p L {pl > 
0), ^r, Pr (pr < 0), u, and \i f° r the constraints ([SD, ©, ©, and (|S]), respectively, for small 
V the minimization amounts to solve a linear system of equations, which possesses a solution 
A = 1 + 0(V 2 ), A n = 0(V) for n ^ 0, A L , P = 0(V), X R , P = 0(V), a = E + 0(V), and 

Xi = o(v). 

The quantites entering the minimization problem within the linear response approxima- 
tions are 



T n (q,p) = £ (* n | a\^a q+r , a |* ) e^ 2 ^, (12) 

r,cr 

^n(0 = <*nh'j |*0) , (13) 

W n = (* n |W|* Q >. (14) 



For the ground state (n = 0), JF (g,p) = fo(q,p). 

We shall assume (a fact justified for the models considered here) that in the absence of 
perturbation the Hamiltonian of the system H is invariant under inversion, and the ground 
state \l/o is even and carries no current, Jo(l) = 0. Owing to the spatial inversion of H, the 
eigenstates \l/ n are either even (\I> 9 ) or odd fif u ). In view of the space inversion, is it natural to 
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choose boundaries located symmetrically (qt = — qa), a mesh comprising symmetric values 
of positive and negative values of p (pl = —Pr = p > 0), and an antisymmetric applied 
potential profile (</?_; = —ipi). The quantities (TT2|) . (IT51) . and (fill) possess the following 
symmetry properties: 

Fg (to, P) = Fg (iR, ~P) , ?u (ql, P) = -Tu («h, -P) , (15) 

W) = -Jg(-l), Ml) = M-l), (16) 
W g = 0. (17) 

Straightforward algebra shows that the equations for even (g) and odd (u) eigenstates de- 
couple, and only odd (u) eigenstates contribute to the current 

J = Ji = 2^2lmJ u (l)lmA u . (18) 

u 

The minimization yields the following linear equations for the Lagrange multipliers X(p) = 
[Xl(p) - Xr(-p)] /2 and X (j) = [xU) ~ xH)} /2 (j,j',P,p' > 0) 

E x (pO E ~p ~~pT (at, pO ^ p) + K fe, p') ^« (to, p)] 

+ E X(/)E^ ^[^(j")^(to,p)+X:(/)^(g L ,p)] (19) 

1 w« 
= o E p — V ^ (to, p) + (&, p)] • 



£>W£ ^M^fa^ + M f x(f) E i^b)'^ (/) = 0[ (20) 



Once they are determined, the relevant expansion coefficients can be computed as 
A,, 



w M ~ 1 

" ■ E A (p)-^(to,p)+ E x(j)^(j) 



2 /< ./• I 



(21) 



which enables to determine the electric current via Eq. ([TBI . Notice that because W u enters 
linearly Eqs. (|T9|) . (1201 . and (l2~Tj) . the current computed via Eq. (fl8|) is proportional to the 
applied bias V. That is, the minimization procedure described above indeed yields a solution 
corresponding to the linear response limit. 



III. MODEL 



The method exposed in Sec. [TT] is general. It will be applied to a concrete system. The 
physical system considered in this paper consists of a single quantum dot connected to two 
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noninteracting leads. It can be described by the Anderson impurity model 

-M L -M L +1 
H = £ lY,J2 a l,* a l,a ~tL J2 Y ( a l,<r a l-l,a + 4-l,cr%cr) 

l=-l CT i = _l (7 

Mr Mr-1 

+ £ r Y Y a L a i,a - tR Y Y { a L a i+i,a + a}+i,cr a i,cr) 

1 = 1 cr 1 = 1 a 

(22) 

~td,L Y ( -l,A + d U-l,a) ~ td,R Y ( a +l,<A + 4 a +l,cr) 

cr cr 

+S g J2 dUa + Ufl d,^d,l, 

cr 

where a^ a {a\ a ) denote creation (annihilation) operators for electrons of spin a in the leads, 
d a = ao,cr {d\ = ao CT ) creates (destroys) electrons in the QD, and n^a = d\d a are electron 
occupancies per spin direction. We shall consider M L = M R = M, t L = t R = t, e L = Er 
(chosen as zero energy hereafter), t djL = t^R = t d , and this ensures the spatial inversion 
assumed above. The dot energy e g can be tuned by varying a gate potential. In view of the 
particle-hole symmetry of model (|22|) . (e g = —U/2 is the particle-hole symmetric point), we 
can restrict ourselves to the range e g > —U/2. The number of electrons iV will be assumed 
equal to the number of sites (N = 2M + 1). The electric potential entering Eq. (flOj) will be 
assumed constant within the reservoirs tpi = —<p-i = — V/2 for < I < M and zero at the 
dot location, ip = 0. 

In isolated electrodes (td = 0), the single-particle energies lie symmetrically around e = 0. 
Therefore, in order to eliminate energy corrections ~ t/M in small clusters, it is advantageous 
to consider reservoirs with an odd number of sites M.— ^ In the noninteracting case (U = 0), 
the resonant tunneling at the Fermi energy (fi = 0) is favored for odd M. In the presence of 
interaction, the Kondo effect occurs when the dot spin couples with electrons of the leads 
at the Fermi level. This coupling is favored if the leads possess single electron levels of zero 
energy, and this is the case for reservoirs consisting of an odd number (M) of sites. 

Although the numerical values used for the results presented in this or in the next Sec. 
IIVI do not represent a special choice, we prefer to employ values corresponding to a possible 
physical realization of the above model, namely chains of QDs of silver. Such QDs were 
experimentally fabricated.— i2L2&22d£iiiid2iM Their properties can be tuned in wide ranges by 
varying the dot diameter (2R) and interdot spacing (D) ^^L&L&LdQ&LdZ.M.^ and the parame- 
ters are well documented in a series of experimental and theoretical works.— liSiMii^ For inter- 
dot spacings up to say, D/2R ^ 1.3, electron correlations are not important^ 1 ^ 1 ^. There- 
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fore, the reservoirs can be modeled, e. g., by chains of nearly touching QDs (D/2R ~ 1.07, 
t = leV). 

In principle, the boundaries could be chosen at < qj? = —q^ < M. Unless otherwise 
specified, we shall choose qh,R = ^2. It amounts to consider the central part to consist 
of the QD and two additional sites, one at either side of the "device", represented by the 
QD of interest. This bears the most resemblance to usual quantum chemical approaches to 
molecular transport, wherein the smallest possible parts of electrodes are included in the 
central part. In Eq. ([2]), there are n QL = 2{M + q L ) + 1 values of r. such that the values of 
site indices are qL ± r = 0, ±1, . . . , ±(M + q^). According to the careful analysis of Ref. [30, 



the values of p in Eq. (T5]) should span the first Brillouin zone of the reciprocal space of 2r. 



IV. CONDUCTANCE THROUGH A POINT CONTACT IN THE ABSENCE OF 
CORRELATIONS 

Let us start with the noninteracting case, amounting to switch off the Coulomb interaction 
(U = 0) in Eq. (122]) . This represents the textbook case of conduction through a single level 
system.— By approaching the resonance, e g — > 0, the transmission becomes perfect. The 
curve of the conductance G(e g ) exhibits a peak characterized by a height G(0) = Go and a 
half-width parameter T = 2t 2 d /t. 

For the numerical calculations in this section, we shall choose a value td/t = 0.4. This 
corresponds to an Ag-QD chain, with a QD in the middle slightly more distant (D'/2R = 
1.24) from its two neighbors than the other QDs in the chains (D/2R = 1.07, cf. Sec. IHIj) . 

In the case of identical reservoirs and contacts, the zero-bias conductance can be obtained 
from the Friedel-Langreth sum rule^ 1 ^ 1 ^ 1 ^ 

G/G = sin 2 (7m d /2), (23) 

where G = 2e 2 /h is the conductance quantum. The dot occupancy per spin direction 
n d = Do-^ol^o-l^o) wn l De computed by numerical exact diagonalization. 

In Fig. [U we present numerical results obtained for 63 sites by means of Eq. ( 1231) . which 
show a peak in the conductance G(e g ) with a half width at half maximum in very good 
agreement with the formula T = 2t 2 d /t. These results are in accord with the fact that 
the single-particle quantum tunneling constitutes the underlying physical phenomenon: At 
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sufficiently large values of e g , the curve for \ogG(e g ) in Fig. [TJ varies linearly with £^ 2 , as 
expected for the transmission coefficient through an energy barrier e g . The lowest excitation 
energy, also shown in Fig. [IJ displays a similar dependence, which confirms that it plays the 
role of a tunneling splitting energy. 

The curve of conductance for chains with 63 sites depicted in Fig. [1] is very close to the 
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. It has been shown there that 



exact result for infinite chains presented in Fig. 5 of Ref. 
the latter result agrees very well with the time-dependent-DMRG result for 64 sites.— The 
size N = 63 chosen by us is the closest to chain size N = 64 used in the time- dependent 
DMRG-calculations compatible with with our choice (N = 2M + 1, with odd M, cf. Sec. 

urn. 




0.01 1 

0.2 0.4 0.6 0.8 e i/2 1 

s 



FIG. 1: (Color online) Results for the normalized conductance G/Gq and lowest excitation energy 
Set (tunneling splitting) obtained by exact diagonalization in the uncorrelated system (U = 0) for 
chains with 63 sites, t = 1, and = 0.4. Note the logarithmic scale on the ordinate and the square 
root of the dot energy(=energy barrier) e g on the abscisa. The linearity of the two curves at larger 
e g confirms the interpretation within the tunneling model. 

We shall now present the results of the SWF-method discussed in Sec. [Til From the point 
of view of the computation time, large chain sizes are more prohibitive for the SWF-method 
than for exact diagonalization. Indeed, by inspecting Eqs. (fl8|) . (Tl9|) . (EDI) , and (12T1) one 
can see that the computing time for J scales as N 6 in the noninteracting case, where the 
number of relevant excitations scales as N 2 . Because the Hamiltonian (|22l) is quadratic in 
the noninteracting case, exact numerical diagonalization can be straightforwardly carried 
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out even for very long chains, e. g., much longer than those that can be handled by the 
time-dependent DMRG. 

Fig. [2] represents a counterpart of the bottom panel of Fig. 1 of Ref . 
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and illustrates a 



characteristic feature of the SWF-method discussed in Sec. [TTJ the reservoirs only constrain 
the Wigner distribution function of incoming electrons, while the distribution of outgoing 
electrons is free. For instance, at the right interface 5f(qR,p) = f(qR,p) — fo(QR,p) is zero 
for p < (in accord with Eq. (jSJ)) but has nonvanishing values for p > 0. Curves for the 
latter case are shown in Fig. [2j depicted for several positive values of p. 
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FIG. 2: (Color online) The difference Sf(qR,p) = f(qR,p) — fo(qR,p) (in arbitrary units) plotted 
versus gate potential e g for 63 sites, t = 1, t& = 0.4. The nonvanishing values at the right interface 
qR and positive momenta p indicate that the Wigner distribution function of outgoing electrons is 
free. The values of k (p = kn /n qL ) are given in the legend. 

We shall now compare the exact results with those of the SWF-method. In Fig. [31 the 
SFW-curve for zero-bias conductance is plotted along with the exact curve. As one can 
clearly see there, the SFW-curve looks completely different, bearing no resemblance with 
the exact curve. Most unphysically, the SWF-conductance vanishes for resonant tunneling 
(e g = 0), where it should attain the maximum value G = Go- 
To compute the SWF-curve of Fig. [31 we have chosen = —2 and qR = 2. A nontrivial 
realistic ab initio calculation is so demanding, that, besides the device (a molecule or a few 
QDs), if at all, only small parts of the electrodes can be accounted for in the heaviest part 
of the computation. Therefore, practically no or very limited freedom remains to choose 
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FIG. 3: (Color online) Results on the normalized conductance G/Gq obtained by exact diagonal- 
ization and the method based on the Wigner function for 63 sites in the absence of correlations, 
U = (same parameters as in Fig. [1]). 

the boundaries qi,R- The fact that in the present case the eigenvalue problem can be solved 
exactly for large systems S enables us to flexibly change and to inspect the impact on 
the solution, and, maybe to make it having some resemblance to the exact one. In Fig. HI 
we present results derived by choosing different qi,R- For all choices, the optimization yields 
minimum values of the total energy £, Eq. (OQ), below the value £q corresponding to the 
trivial situation A n = 0, Vn ^ 0, i. e., the "condensation" energy 5W = £ — £ is negative. 
However, the conductance changes only by factors of the order of unity. Definitely, it cannot 
be made more akin to the exact G. 

The curves presented in the above Figs. El El and H] have been deduced by constraining 
the solution to satisfy the equation of continuity, as discussed in Sec. [Ill In Fig. EJ results 
obtained by imposing the equation of continuity (depicted by the line denoted by uniform) 
are compared with those derived without imposing this equation. As visible there, the values 
of the electric current do display a significant dependence on site. The latter is indicated 
by the numbers in the legend of Fig. [5j Neither the value of the curent through the QD 
nor the average along the chain (label and average, respectively) coincides or reasonably 
approximates that deduced by imposing the equation of continuity. It was claimed that, 
although in principle necessary because of certain approximations (see Sec. IVII for more 
details), there was no need to impose the equation of continuity for the calculations reported 
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FIG. 4: (Color online) Results for the conductance G/Gq and the bias-induced condensation energy 
5W (the latter in arbitrary units) obtained by means of the SWF-method for boundaries chosen 
a t Qr = —Ql = 2, 3, 4, 5 (values increasing upwards for the G-curves and downwards for the 5W- 
curves at, say, e g = 0.5). Noteworthy, the conductance does not sensitively depend on the choice 
of boundaries. 



in Refs. 
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Obviously, the fact that in our case the minimization without imposing 
Eq. (jHJ) leads to a solution for \1> that violates the equation of continuity is not the result 
of any approximation: except for the SWF-method itself, our results are affected by no 
further approximation. To conclude, Fig. demonstrates that the SWF-method does not 
automatically satisfy the continuity equation, not even approximately. 



V. CONDUCTANCE THROUGH A POINT CONTACT IN THE PRESENCE OF 
CORRELATIONS 

In this section we shall apply the SWF-method exposed in Sec. [Ill for the case of nonvan- 
ishing U, where correlations are known to play an important role. The physics of model ( !22l 
with U 7^ is also well understood.— For strong interaction (U) and low temperatures, by 
varying the gate potential e g , one observes plateaus of well defined dot charge, corresponding 
to a dot that is empty, singly, and doubly occupied: rid = (0 < Ed), = 1 (-U < Sd < 0), 
and rid — 2 (e^ < —U), respectively. This behavior can be demonstrated by exact numerical 
diagonalization in small clusters, as illustrated by the curves of Fig. [61 
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FIG. 5: (Color online) SFW-results for the conductance G/Gq of the chain with 63 sites at t = 1 
and td = 0.4. The curves computed without imposing the continuity equation exhibit a strong site- 
dependent current and substantially depart from that labelled by uniform, computed by imposing 
this equation. Nor the average taken along the chain (label average) represents a satisfactory 
approximation of the latter. The number of site q in the chain is specified in the legend. The dot 
is located at q = 0. 

For low temperatures but above the so-called Kondo temperature Tk (often much smaller 
than IK), the conductance G(e g ) exhibits two narrow Coulomb blockade peaks located at 
e g = —U and Eg = 0, while in between it almostly vanishes ("Coulomb valley"). By 
decreasing the temperature below Tk, correlation effects in the singly occupied state yields 
a sharp (Kondo) resonance in the density of states at the Fermi level, and this gives rise 
to a characteristic plateau of width ~ U in the curve of G versus e g . In the middle of the 
Coulomb valley (e g = —U/2), perfect transmission occurs, leading to the ideal conductance 
value G = 2e 2 /h (unitary limit). 



Numerous results obtained for the model f l22l) 



published in the literature; see, e. g., Refs. 



35 



jy considering semi-infinite leads were 
. A comparison between the results for 



semi-infinite leads and the SWF-method would make no sense. At zero temperature, the 
case for which the SWF-method was developed, in the range —U ^ e g < for the realistic 
case of semi-infinite electrodes the conductance is dominated by the Kondo plateau. Its 
formation requires chain sizes larger than the Kondo cloud, which extends over a number of 
sites £k ~ t/Tjc- The latter rapidly grows (exponentially for large U) beyond the sizes, which 
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neither the exact diagonalization, nor the SWF, or often even the DMRG approach^ can 
handle. No Kondo peak can be formed for chains shorter than the Kondo screening length 
£ K . However, in short chains the G(e 9 )-curve should still display the Coulomb blockade 
peaks at e g = and e g = — U. 

Although the clusters considered in this section comprise a small number of sites, the 
results are significant. Since our main purpose here is to address the issue of the validity of 
the SWF-method, we do not intend to discuss finite-size effects here. However, we do not 
expect that they are essential: the inspection of Fig. [6] reveals that the differences between 
chains with seven (M = 3) and eleven (M = 5) sites are not substantial. 

The primary reason why we restrict ourselves to chains with seven sites is technical, but 
this also rises supplementary doubts on the applicability of the SWF-methods for systems 
of interest for molecular electronics. Eq. (ITT1) is an expansion over the complete set of eigen- 
states of H and, according to Eq. (Tl8l) . in principle all eigenstates of odd parity contribute 
to the current. Calculations contradict the naive expectation that in the linear response 
approximation only a reduced number of excited states are important. For the couplings 
U employed in Fig. [3, out of a total of 1225 states with spin projection S z = +1/2, the 
first 300 states are not always enough to reach convergence. For eleven-site chains (the next 
larger size of interest), there are 213444 eigenstates with S z = +1/2. One would probably 
need to target many thousands thereof in order to get the matrix elements necessary for 
convergent results. For this formidable task, one should run the Lanczos procedure three 
times, and this separately for each of the matrix elements entering Eqs. ( fl2|) . ( 1T31) . and ( fl4T) . 
in a manner similar to but more involved than that employed to compute frequencies and 
intensities of optical lines.— illi^i^ The method of only computing convoluted spectra by 
means of the continued fraction algorithm^^^ is inadequate for this purpose. 

Our results on the conductance computed by means of the SWF-method are collected in 
Fig. (Notice that only the halves of the curves situated at the right of the particle-hole 
symmetric point e g + U/2 = are shown.) As one can see there, the curves for the con- 
ductance, calculated for several values of U and td, exhibit maxima at the gate potential 
values e g > (e g + U/2 > U/2), i. e., at the position where the Coulomb blockade peaks 
are expected. Their distance from the ideal Coulomb blockade location increases with t^, a 
behavior similar to that of their width. While this behavior is physically plausible, unfor- 
tunately, the prediction of the SWF-method for the height of the peaks is quite unphysical: 
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FIG. 6: (Color online) Dot occupancy computed by exact diagonalization for 7- and 11-site clus- 
ters (solid lines and points, respectively) plotted versus dot energy e g . Notice the small difference 
between the results for 7 and 11 sites, indicating that finite-size effects play a reduced role. 

the height is found to vary roughly inversely proportional to td- The consequence of this de- 
pendence is that, as visible in Fig. [3, the height of the G-peak even attains values exceeding 
the ideal value G . To conclude this section, the results of the SWF-method are unphysical; 
the conductance of correlated systems evaluated by this method cannot be trusted. 

VI. DISCUSSION 

We have two comments on the SFW-method. They concern this method in general, and 
are not related to the linear response limit. Firstly, in Ref. |2l|, it was claimed that the impo- 
sition of Eqs. (jHJ) is necessary, because the continuity equation dp/dt + \7-j = 0, derived from 
the Schrodinger equation in the presence of local interactions, does not hold for nonlocal 
interactions. Accordingly, position-dependent currents could be an effect of ab initio molec- 
ular electronics calculations employing nonlocal effective core potentials or pseudopotentials, 
or a result of truncating the molecular orbital basis set or the CI (configuration interaction). 
While all these may in general be sources of violating the continuity equation, for the SWF- 
method there still exists another reason to impose the constrains (jSJ) in a steady state. The 
wave function \l/ determined by means of the SWF-method does not represent an eigenstate 
of the Hamiltonian. \l/ is time independent only because the formalism is time independent, 
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FIG. 7: (Color online) Coulomb blockade peaks of the zero-bias conductance for 7-site clusters 
as predicted by the SWF-method. The solid, dotted, and dashed lines correspond to the values 
td = 0.125, 0.25, and 0.5 respectively. The values of U are given in the legend, and t = 1. Notice 
that, unphysically, the maximum conductance is predicted to increase with decreasing td even 
beyond the ideal value Go (G/Gq > 1)- 

and not the result e. g., of taking the limit t —>■ oo to get a steady state. No demonstration 
has been given in the works dealing with the SWF-methocPQ^ 1 ^ 1 ^ 1 ^ that the minimiza- 
tion procedure without imposing Eqs. (|Hj) yields a solution compatible with the continuity 
equation. While the imposition of Eqs. (jHJ) is seemingly unnecessary for the case considered 
in Ref. |2l|, there is no rationale for this in general. For illustration, we have presented a 
counter-example in Sec. [TV] The equation of continuity dni/dt = (i/h)[ni, H] = —ji + ji-i 
(cf. Ref. I33I) holds for the model Hamiltonian (j22j) . and nevertheless the current computed 
without imposing Eqs. (jHJ) is strongly site-dependent. 

The second and more important point concerns the manner of imposing boundary condi- 



tions in Ref. 
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J. In the approaches based on the Wigner function within the single-particle 
approximation the influence of the applied electric field is accounted for both at the boundary 
conditions, via the shift fii — = (cf. Eqs. (JBJ) and PJ), and on the electron dynamics 
within the device.— Within the methods based on the NEGF the applied voltage is usu- 
ally considered solely via the the chemical potential imbalance, and many-body^ methods 
are employed to treat correlation effects due to interactions within the device without ap- 
plied field. In both cases, the applied bias represents the driving force of current flow. In the 
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SWF-method, the effect of the applied field at the boundaries is entirely neglected (cf. Eqs. 
(jSJ) and (jSJ)). Let us suppose that (i) we would be able to reliably solve the minimization 
problem as prescribed by the SWF-method and exactly (or at least very accurately) deter- 
mine \I> for very large systems, including large parts of the reservoirs (assumed identical), 
and (ii) in the latter the single-particle description applies. Then, according to Eqs. (jSJ) and 
dHl), the Wigner functions at the boundaries would reduce to the Fermi distribution functions 
of the left and right reservoirs characterized by the same chemical potential \i. This would 
imply that there would be no difference between incoming electrons from the left and right 
reservoirs. Then, it is not at all surprising e. g. that in the extreme case, where all sites are 
noninteracting and identical (U = and = t), instead of being maximum, G = Go, the 
conductance vanishes for e g = 0, as visible for the SWF-curves of Figs. [21 HI andEJ In reality, 
the correct wave function ^ describing the steady state current flow should yield a Wigner 
function that reduces at the boundaries to the Fermi distribution functions characterized by 
different chemical potentials ^l,r = /i ± eV/2. 

Since the above analysis reveals that the boundary conditions (j5J) and (jEJ) are inadequate, 
attempting to mend the SWF-method would be desirable. In view of the above consid- 
erations, perhaps the most natural attempt would be to modify Eqs. (jSJ) and dm by using 
instead of \l/o the ground state $o of the system in the presence of the applied potential, i. e., 
H T \&o) = E T)0 \& ). Although the Wigner functions entering the r.h.s. of Eqs. (jSJ) and 
calculated by using $o istead of \l/o, do not necesarily reduce to the left and right Fermi dis- 
tributions, this procedure would at least account for the chemical potential imbalance at the 
boundaries. However, as revealed by a straightforward analysis, this modification does not 
yield the desired improvement: the solution of the minimization is just \I/ = $o- This solution 
obviously satisfies the boundary constraints (jSj) and dm as well as the continuity equation 
(electrical perturbations only depend on electron density), and, being the ground state of 
Ht, it trivially minimizes £ of Eq. [TJ This result is obviously general, i. e., it holds beyond 
the linear response approximation. Still, as a verification, we have performed the modifica- 
tion \l/o — > in the r.h.s. of Eqs. (jSJ) and (ED, and carried out straightforward calculations 
within the linear response approximation. They yield \^/) = \^f ) +J2 n ^o W„/ (E — E n )\ty n ), 
and in the r.h.s. one immediately recognizes the ground state <3>o of Ht = H + W in the 
first-order of perturbation theory. Unless the system is superconducting, the current (con- 
ductance) vanishes in the state described by the wave function ^ = $ . So, even with this 
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"remedy" , the SWF-method is unable to describe electric transport through nanosystems. 



VII. CONCLUSION 

Several recent studies proposed and applied a many-body time-independent method 
to compute steady-state electric transport in molecular systems, whose key ingredient 
was the formulation of the boundary conditions in terms of the Wigner distribution 
function.— The fact that this approach yielded values of the electric current through 
molecules comparable with those measured in experiments, which are usually orders of mag- 
nitudes lower than the predictions of other theoretical treatments, was considered very 
encouraging. However, the mere fact that a theoretical method compares favorably with ex- 
periments cannot be taken as support for its correctness. It should also be able to correctly 
reproduce well established results. 

In this paper, we have presented results demonstrating that the SWF-method is unable 
to reliably evaluate the zero-bias conductance of the simplest uncorrelated and correlated 
systems of interest for molecular and nanoscopic systems, namely a single QD. It fails to 
retrieve the result G = Go for resonant tunneling through a single QD without correlations, 
where it predicts a vanishing conductance instead. In the presence of correlations, the 
conductance at the peaks of Coulomb blockade is unphysically predicted to increase with 
decreasing dot-electrode coupling (to) and can even exceed the conductance quantum Go. 

While the idea of formulating boundary conditions in terms of the Wigner function for 
correlated many-body systems is interesting, the manner in which it was imposed in Ref. 
20l turns out to be inappropriate. It misses the fact that, in accord with our physical 



understanding, the current flow is due to an asymmetric injection of electrons from reservoirs 
into the device, and that injected electrons are very well described by Fermi distributions 
with different chemical potentials. Moreover, as results from the analysis at the end of Sec. 
IVIt unfortunately there is no simple remedy of the SWF-method; the modification of the 
boundary conditions in the spirit of Ref. such as to account for a nonvanishing chemical 
potential shift does not yield the desired improvement. 

In addition, as a side note, we believe that the very fact that astonishingly numerous 
states with high excitation energies are found to contribute to the conductance in the linear 
approximation is an indication that the SWF-method, even if it were physically sound, 
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would be of little pragmatical use for strongly correlated nanosystems. Because it would be 
hardly conceivable that ab initio calculations for real molecular systems, by far more complex 
than the presently considered model, could provide a wave function \1/ with the accuracy 
needed for reaching reliable convergent results. These considerations raise doubts on the 
possibility to develop a viable SWF-approach for the electric transport through nanoscopic 
and molecular systems. 
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